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Abstract We propose an adaptive control law that allows one to identify unstable steady 
states of the open-loop system in the single-species chemostat model without the knowledge 
of the growth function. We then show how one can use this control law to trace out (recon- 
struct) the whole graph of the growth function. The process of tracing out the graph can be 
performed either continuously or step-wise. We present and compare both approaches. Even 
in the case of two species in competition, which is not directly accessible with our approach 
due to lack of controllability, feedback control improves identifiability of the non-dominant 
growth rate. 

Keywords Chemostat, growth, identification, competition, slow-fast systems, numerical 
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1 Introduction 

We recall the classical chemostat model [31] for a single species (biomass b) consuming a 
substrate (mass s): 
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s = -fi(s)b + D(s in - s 
b = n\s)b-Db 



(1) 
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where the dilution rate D (the input) is the manipulated variable, which takes values in 
a bounded positive interval [Z5 m ; n ,D max ], and /i(-) is a non-negative Lipschitz continuous 
function with /i(0) = . 

We consider here the following scenario: the function /i(-) is unknown and possibly 
non-monotonic. Our objective is to reconstruct the graph of the function fi(-) on the domain 
(0,jj n ) by varying the input D in time. On-line measurements are only available for the 
variable s (that is, s is the output). This setup is realistic for experimental investigations such 
as in [3], however, demonstrations in this paper are based entirely on simulations of models 
such as system (1). The present paper analyzes and expands the ideas initially proposed by 
the authors in the conference paper [30]. 

Remark: Using model (1) tacitly assumes that the yield coefficient of the bio-conversion 
is known. This is why [i(s)b appears with the same pre-factor 1 (once positive, and once 
negative) in both equations of ( 1 ) without loss of generality. 

The problem of kinetics estimation in biological and biochemical models has been 
widely addressed in the literature ([11,12,13,5,14,23,10]), either as a parameter estima- 
tion problem (one chooses a priori an analytical expression of the function /i(-)), or as an 
on-line estimation of the kinetics (one aims at determining jl(s(t)) at the current time t). 
The theoretical identifiability of the graph of fi(-) has been thoroughly studied in [8]. In this 
paper, a practical method has been proposed to reconstruct the graph of based on a 
Kalman observer under the approximation that the function z(f) = fi(s(t))b(t) has a third 
time derivative equal to zero. 

Here, we propose a different method that does not make any approximation of the dy- 
namics. Our method exploits that it is sufficient to find the complete branch of equilibria of 
system (1) to identify the graph /i(-). This reduces the system identification problem to a 
combination of two problems: finding the equilibria (a root-finding problem) and stabiliz- 
ing them (a feedback control problem). Both of the latter two problems are in theory easily 
solvable with standard methods as we will explain and illustrate in Sections 2-5. The most 
difficult obstacle in practice is the implementation of a real-time feedback loop measuring s 
and adapting D with sufficient accuracy. 

When the growth function is monotonic, a common way to reconstruct points on the 
graph of the growth function jU(-) is to design a series of experiments fixing the dilution 
rate D with different values and wait until the system settles to a steady state (s*,b*) [5]. As 
long as D is less than [i(s m ), it is well known that the dynamics converges to a unique posi- 
tive equilibrium that satisfies fi(s*) = D (see for instance [31]). This technique requires the 
steady state to be stable in open loop, and consequently cannot reconstruct any part of the 
graph of a function where pt is non-increasing (such as the example, shown schemati- 
cally in Figure 1). Furthermore, the global convergence of this method is not satisfied in case 
of bi-stability, which is present in (1), with non-monotonic growth functions pt (see again 
[31]). 

An alternative approach is to fix a value of s, say s and design an adaptive control law 
£>(•) that stabilizes the system about the steady state (s, s m — s), with the value of D con- 
verging to n(s). Several adaptive control laws have been proposed in the literature for this 
problem. Nonlinear feedbacks require the knowledge of the growth function fi(-), such as 
linearizing controls [5,20] Extensions that are robust with respect to uncertainty on /i(-) 
have been proposed [25] but do not provide the precise reconstruction of H(s). Several non- 
linear PI based controllers, that do not require the precise knowledge of fi(-), have been also 
proposed [26,28], but saturation and windup is often an issue (see [18, 19] in similar frame- 
works). In [3], a dynamical output feedback has been proposed to globally stabilize such 
dynamics without the knowledge of /J. ( ■ ) and under the constraint D 6 [Anin, Anax], but it 
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Fig. 1 Domains of stability and instability in open-loop 



requires the growth function to be monotonic. More recently, a saturated PI controller cou- 
pled with an observer, dedicated to the non-monotonic case, have been proposed to stabilize 
the dynamics about a nominal point that maximizes the biomass production [29]. 

In this work, we first propose that it can be useful to introduce a feedback control loop 
into (1) to identify the growth function [1 of the open-loop system (1) (that is, (1) with 
constant input D). The feedback control law is initially a simple saturated proportionate 
controller: 

D(s,D,s) = sat [zwzw] (D - d (a - S)) , (2) 

where D and s are reference values, and G\ > is the linear control gain. To ensure realistic 
values for the input D, feedback law (2) encloses the linear feedback rule into the saturation 
function 

^max if A: > .D m ax> 

sat [D rain ,D raax ] (*) = <* if x e [D min , £> max ] , 

t ^min if X <C D m [ n , 

where the limits D m [ n and Z) max are the extreme dilution rates that can be achieved experi- 
mentally. In sections 3 and 5 we will then explore adaptation rules for the reference values 
(f, D) which ensure that asymptotically for t — > °° the input satisfies 

D(s,D,s)=D (3) 

or, equivalently, the output satisfies 

s = s, (4) 

such that the control is asymptotically non-invasive. The result is a new adaptive control 
law that stabilizes the dynamics about any desired equilibrium point without requiring a- 
priori knowledge of its location, and whatever is the monotonicity of the growth function. 
One requirement on the adaptive law is that it should work uniformly well around a local 
maximum of /i (non-invasive feedback laws such as filtered feedback [ 1 ] or time-delayed 
feedback [24] do not achieve this). 

We note that our adaptation rules will be much simpler than classical adaptive control 
laws [5,4]. Usually, adaptive control aims to achieve a desired output regardless of changes 
in the underlying system. We only adapt the reference values to make the control input 
vanish and find branches of equilibria and bifurcations of the underlying system, similar to 
numerical continuation [2]. While classical adaptive control requires system identification 
(an inverse problem) at some stage, our adaptation solves a root-finding problem, which is 
simpler. 

Throughout our paper we assume that the output s can be sampled and the input D 
can be adjusted in quasi real-time. If the sampling period T is not negligible, the approach 
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Fig. 2 Sketch illustrating shape of growth function n and set of admissible reference values 



presented here can still be applied. However, one then faces the problem that feedback sta- 
bilization of an unstable equilibrium at output s becomes sensitive with amplification factor 
~ exp(— p.'(s)T) to disturbances (note that fi'(s) < for unstable equilibria). 

We show in Sections 4 and 5 that the feedback law (2) can be combined with an adaption 
rule for (s,D) to reconstruct the graph of the growth function, even in the case of non- 
monotonic growth functions. Section 4 presents a dynamic adaption, whereas Section 5 
introduces a step-wise adaptation. In Section 6 we investigate the case of two species that 
compete for the same common substrate. 



2 Global stability of the simple feedback law (2) 

Let us first prove that the feedback law (2) is, within reasonable limits, globally stabilizing. 
Suppose that we choose the reference value s from an interval [.? mm ,Sin) C (0,i; n ), and that 
the limits on the input cover the growth function fi on this interval: 

D min <n(s) for alls 6 [.s 

min j ^in ], (5) 
Anax > MOO for all s e [0, J in ] . (6) 

These conditions mean that the graph of /i does not cross the thick parts of the horizontal 
lines D m i a and D max bounding the grey area in Figure 2 from below and above. 

Proposition 1 Suppose that the reference values (s,D) in feedback law (2), 

are chosen from the rectangle [■s m in,Jin) x \P mim^ max], that the growth function pi satisfies 
(5)-(6), and that the gain G\ is chosen sufficiently large, that is, 

Gi > — min fi (s), and (7) 
se[0, Si „] 



pl{s in )-D 



(8) 



Then the controlled system (1) with D = D(s,D,s) has a stable equilibrium (s e „,b e n) G 
[0,Jj n ) x (0,°°), which attracts all initial conditions (j(0),fc(0)) 6 [0,J; n ) X (0,°°). 
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Proof If D > 0, and the growth function fi satisfies /i(0) =0 and, for s > 0, fi(s) > then 
the set 

R = {(s,b):se[0,s in ),b>0} 

is positively invariant (that is, trajectories starting in R will stay in R for all positive times). 
Furthermore, all trajectories starting in R approach the subspace (called stochiometric set in 
[29]) 

T = {(s,b)eR:s + b = s iD } 

with rate at least Z) m j n forward in time. This implies that it is sufficient to check if all trajec- 
tories in T converge to a unique equilibrium. On T the equation of motion can be expressed 
as a differential equation for s only: 

i=p(sA^-li(s)]K-i]=[satp^^](D-(h(s-s))-ii(s)] [s m -s]. (9) 

First, let us check that the equilibrium at s = s m is unstable. The term —G\ (s m — s) is negative 
such that D — G\ (s- m — s ) < D max for all admissible D. Assumption (8) guarantees that D — 
G\ (sin — s) < n(s- m ). Assumption (5) guarantees that also D m in < M( J in)- Hence, 

D{s in ,D,s)-n(s in ) <0 

for all admissible (s,D). Thus, the prefactor of S[ n — s in (9) is negative such that the equi- 
librium at is unstable for all admissible (s,D). 

Since s > at s = 0, there must be other equilibria of (9) in (0,Sj n ), which are given as 
solutions s eq of D(s eq ,D,s) = /i(s eq ). Now let us check indirectly that none of the equilibria 
can satisfy D min = n(s eq ). 

Assume that (9) had an equilibrium j eq with D m i n = /i(s e q). Then s eq has to be less than 
Smin due to assumption (5). However, if .s eq < s m [ n , then D — G\ (s eq — s) > D > D m in for 
all admissible (s,D). Hence D(s eq ,D,s) > -D m j n (recall that D mm = fJ-{s eq ) by assumption 
of the indirect proof) such that D(s eq ,D,s) — n(s e q) > 0, which means that s eq cannot be 
equilibrium, establishing the contradiction. 

Assumption (6) excludes that equilibria with fi(s eq ) = D max exist, hence all remaining 
equilibria s eq 6 (0,i; n ) must satisfy 

D-Gl(s eq -s) =/i(s e q)- 

Condition (7) ensures that this equation has a unique solution and that this solution corre- 
sponds to a stable equilibrium (which must be in (0,Ji n ) because the boundaries of (0,.Sj n ) 
are inflowing for (9)). □ 

Proposition 1 ensures that the output Seq of the controlled system (1) with (2), after 
transients have decayed, is a well-defined smooth function of the parameters (s, D) as long as 
(s.D) are chosen from (i m i n ,i ma x) X (Anin, Anax)- We express this fact by using the bracket 
notation: 

s eq (s,D) = lim.s(r) where s is output of (1), (2). 

The function j- e q can be evaluated at any admissible point by setting the parameters (s,D) 
in the definition (2) of the feedback rule, waiting until the transients of ( 1 ) have settled, and 
then reading off the output s. Equilibria of the uncontrolled system can then, according to 
(4), be found as roots of s eq (s,D) — s. More specifically, we know that, for any admissible s, 

D = n(s) if and only if s eq (s, D) = s. (10) 

Relation (10) permits us to identify fi(s) as the unique root of s eq (s, •) — s. Sections 3-5 will 
explore two strategies to find this root for a range of admissible s efficiently. 
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3 An adaptive control scheme 

The first strategy is a dynamic feedback that comes on top of the feedback law (2) for D. We 
treat D not as a parameter but introduce an additional dynamical equation for D, achieving 
local convergence of the output s to any reference value s£ (0,s; n ) without the knowledge 
of the growth function /x . Then the asymptotic value of D allows one to reconstruct the value 

m> 

Proposition 2 Fix a number s 6 (0, s m ) and take numbers Z) m j n , D mdx that fulfill < D m \ n < 
fj.(s) < D max . Then the dynamical feedback law 

D(s,D)=sat [Dm<M (D-G^s-s)) 
D=-G 2 (s-s)(D-D min )(D mea -D) 

exponentially stabilizes the system (1) locally about (s,b) = (s,S[ n — s), for any positive 
constants (Gi , G 2 ) such that G\ > —fl'(s). Furthermore one has 

lim D(t) = u(s) 

t->+°o 

Remark. Our adaptive control is in this case similar to a classical PI controller. The 
quantity D is playing the role of the I part, but staying bounded by construction. 

Proof Locally about s = s , the closed loop system is equivalent to the three-dimensional 
dynamical system 

s = -Lt(s)b+{D-Gi(s-s))(s in -s) 
b = p.(s)b- (D-Gi{s-s))b 
£> = -G 2 {s- s) (D - D min ) (D max - D) 

This system admits the unique positive equilibrium E* = ( j, s m — s,fi(s)). 

For simplicity, we write the dynamics in the variables (z, s, D) coordinates, where 7 is defined 

as z = s + b: 

z = (D-G l (s-s))(s m -z) 
s = -H (s) {z-s) + (D-Gi (s - s) ) (s m - s) 

D= -G 2 (S-S)(D- D min ) ( Anax " D) 

The Jacobian matrix at E* in these coordinates is 

\ 

-H(s) -{}i'(s) + Gi)(s in -s) s m -s 
-G 2 (fi{s-)-D min ){D mlix -^(s)) / 

Its eigenvalues are X\ = —/J.(s) < and X 2 , A3 as eigenvalues of the sub-matrix 

/ -(/*' (s) + Gi ) {s in - f) s in -s\ 
\-G 2 (li(s)-D min )(D m . dx -Li(s)) J 

Then, one has 

det(M) = G 2 {ti{s)-D min )(D max -ii{s))(s in -s) 
tr(M) =-{ii'{s) + G l ){s m -s) 

and concludes about the exponential stability of E* when G 2 > and G\ > — /i'(f). Finally, 
one obtains from (11) that D or D converges toward the unknown value /J.(s). □ 
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Fig. 3 The behavior of feedback law (11) for two different values of * in the (s,D)-projection. Parameters: 
G\ = 2, Gi = 2, output disturbance defined in (13). 

Note that the assumptions in Proposition 2 (for example, on the gain G\ ) are weaker than 
those of Proposition 1 as Proposition 2 is only concerned with local stability and a single 
reference value s. 

Figure 3 demonstrates how control law (11) stabilizes an equilibrium with output s for s 
in the increasing (left panel) and decreasing (right panel) part of the growth law /i. For our 
single-species demonstration we choose the non-monotonic Haldane function 

*«=i+^Ia? (12) 

and s; n = 1. To illustrate the effect of disturbances, we super-impose a rapid oscillation onto 
the measurements of output s, such that the output has the form 

■^output = s[l + <5cos(3f)sin(?)], where 5 = 0.05. (13) 

The grey background curve in Figure 3 shows n(-), which is clearly non-monotonic on the 
domain (0,Ji n ) (s- m = 1). 



4 Reconstruction of the entire growth function 

Now, we can trace out the whole graph of /i dynamically by letting s change slowly with 
time as solution of the simple dynamics 

s = es(sin-s) (14) 

to explore the right part of the graph of fi(-) when e is a small non-negative number, and to 
explore the left one when e is a small non-positive number. During the reconstruction of the 
graph, the gain Gi has to be been chosen uniformly large according to (7). Figure 4 shows 
how the adaptation rule (14) together with (11) reconstructs the entire graph of the growth 
function. 
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Fig. 4 Dynamical adaptation using (14) to explore dynamically the left and right part of the graph of the 
unknown function n(-). Parameters: for the left panel, e = +0.01, for the right panel, e = —0.01, otherwise 
identical to Fig. 3. 



5 Step-wise adaptation of the reference values 



In this section, we propose an alternative to the continuous adaptation of D and s: we treat 
the root problem = s eq (s,D) — s with ordinary numerical root-finders such as the New- 
ton iteration. We present here an approach that combines the two steps of the method (the 
adaptive control and the continuation) in a step- wise framework. 

In an experimental setting one will have to adapt the numerical methods to the lower 
accuracy of experimental outputs (see [27] for a demonstration in a mechanical experiment) 
but for this paper we restrict ourselves to a numerical demonstration. In the single-species 
chemostat one profits from the knowledge of an approximate derivative of s eq with respect 
to D, making the Newton iteration more efficient. Suppose, we plan to identify the growth 
function n in a sequence of points = Sq + k8 (where 8 > is small). The function values 
H(Sk) are the roots of Seq(*jti*) — *it- Then we can approximate 



1 



Gl+/i'(*eq(Ast)) 



1 



Gi + 2=^=1 



to obtain the iteration 



starting from D = 1 , or (for k > 2) 



D-Dk-l 



(15) 



= .D*-i + 



Dk-i-Dk-2 



For the initial step (k= 1) the derivative of s eq has to be either guessed or approximated with 
a finite difference (we used the latter in our numerical simulations). 

Note that at no point it is necessary to set the internal states s or b of system (1). Only 
the reference values (s,D) have to be set. Figure 5, panels (a) and (b) shows the output 
of a simulation with the step- wise adaptation using Newton iteration (15). The black dots 
correspond to values at which the control was accepted as non-invasive. By gradually tracing 
out the graph of fi, one achieves small and rapidly decaying transients in every evaluation 
of s eq (which involves running system (1) with control until transients have settled). This is 
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Fig. 5 Simulation of step-wise adaptation. Panels (a) and (b) use (15) where = 0.05k. Panels (c) and (d) 
use ( 16)— ( 17) where S = 0.05, both cases subject to output disturbance (13). Panels (a) and (c) shows the time 
profile of output .s and input D throughout the run. Panels (b) and (d) show the evolution in the (s,D) -plane 
in grey. Black dots indicate when convergence was reached (|.v eq — Sk\ < tol for panels (a) and (b), settling 
of transients for panels (c) and (d)), and the iteration moved on to the next s^. Note the shorter time scales 
on the A-axis of panel (c) compared to panel (a). Parameters: tol = 10~ 3 , Gi = 1, D max = 0.2, D m i„ = 0.02. 
A check if the transients have settled was performed every 20 time units. Transients were accepted as settled 
if the standard deviation of s (std(.v)) on the last interval is no longer smaller than 0.9std(.s) on the previous 
interval. The average over the last interval is used as the resulting equilibrium. 

so because the transients all lie inside the subspace {{s,b) : s + b = s m } after system (1) has 
run at least once. Second, the initial offset from the equilibrium is always small, because the 
adjustments of s and D are small. 

5.1 A simplified step-wise scheme 

The scheme (15) permits one to find n(Sk) for an a-priori prescribed set of admissible ab- 
scissae Sk. If one wants to recover only the graph of pi one does not need to prescribe the 
sequence §k a-priori, thus, avoiding a Newton iteration. Suppose that we know already two 
points Pk-\ = (jjfc-i>Z)fc_i) and pk = {sk,t>k) on the curve (s,n(s)). Then we set 

(Ved.<t+ 1 , D pre d.k+ 1 J = Pk + O r 77 , (16) 

llp*-p*-lll 

where 8 > is the approximate desired distance between points along the curve (s,fl(s)), 
and run the controlled experiment with the reference values (s,D) = (ipi-ed,/t+i,^pred,jfc+i) in 
(2) until the transients have settled to obtain the next point on the curve 

Sk+\ — s eq(Spred,k+l,Dp [e d t k+l) 

Dk+l = 0(4+l>Z)pred,*+l!*pred,*!+l) (17) 
= -Dpredjfc+l — G\ (Sk+\ — .Jpredjt+1 ) 

This simplified procedure cannot guarantee the identification of /i at prescribed equidistantly 
spaced values of s but finds fl(sk) for a (nearly evenly spaced) sequence lj given by the 
intersections of the lines D = D pK ^k — Gi(s — J pre d,yt) with the graph D = n(s). 
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Figure 5, panels (c) and (d), demonstrate the speed-up using the simplified scheme (16)- 
(17) (note the times at the abscissae). The difference to Figure 5(a,b) is that the values s^ at 
which the growth function is evaluated are not exactly equidistantly spaced. 

6 The two species case 

Let us now consider an extension of the chemostat model (1) that considers two species 
which compete for the same substrate. The two-species model can be written as follows 

2 

« = -£iUi( s ) fe i + D ( s in-*) (lg) 

bi = lk{s)bi-Dbi (i=l,2) 

The two-species model has co-existing equilibria E* , which correspond to the state where 
species i is present and the other species 3 — z is suppressed. The following proposition shows 
first that feedback stabilization based on input D and output s breaks down in general for 
the equilibrium corresponding to the species with the smaller growth rate (the suppressed, 
or non-dominant, species). Then we state what eigenvalues the linearizations at equilibria 
have for our specific control laws, (2) and (11). 

Proposition 3 Fix s£ (0, s; n ) and consider the equilibrium E\ — (s, 0, Si n — s). 

1. (Suppressed equilibrium not stabilizable) Let }i\ (s) > }li{s), and D() be a feedback 
D(-) of the form 

D = /(*,§), $=g[s,$), (|eK*J (19) 

with f(s,0) = and g(s,0) = 0. Then the equilibrium Et of system (18) with feed- 
back D is unstable. 

2. (Eigenvalues for simple feedback) Using feedback law (2), 

D{s,D,s) = sat [Dmm:Dm ^ (S-Gi(s-s)) 

with D m i n and Z) max such that D m [ n < j±2 (s) < D max , the linearization of system (I) in the 
equilibrium E\ has the eigenvalues —\Li{s), \l\ {$) — jttzC?) an d — {^(s) + G\)(si n — s). 

3. (Dominant equilibrium stabilized by dynamic feedback) If Hi(s) < fJa(s), then the 
feedback (11) exponentially stabilizes the system (18) locally about Et, for any positive 
constants (Gi ,(72) such that G\ > — /i^(s). Furthermore one has 

^lirnS(f) = n 2 (s) 

Proof Consider the dynamics of the two-species model (18) with feedback D given by (19) 

{S = - £ Mi (*)£>> + /(*. I ) ( J in - *) 
bt = IM(s)bi-f(s,Z)bi (i=l,2) 

We write this system in (z,bi,b2,^) coordinates with z = s + b\ -\-b%: 

i = /(z-*i-*2jI)(%-z) 
k = (LH(z-bi-b 2 )-f(z-b 1 -b 2 ,Z))b i 

i =g(z-bi-bt,S). 
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Point 1: at equilibrium E\, the Jacobian matrix J\ possesses the following form in 

(z, b\ , Z»2, |) coordinates 



(-IH{s) 0\ 
/ii(j)-ft(i)0 
★ * * * 

V * * * -k J 



(20) 



which has the positive eigenvalue p.\ (s) — fJ-2{s)- This proves that £| is unstable whatever 
the choice of the feedback D(-). 

Point 2: we can be more specific about the form of J\ for the simple feedback law (2). Since 
Am n < < Anax, the feedback is in its linear regime, such that d s D = —G\. Thus, 

(component E, is absent) 

/-Ma(i) \ 

y 2 *= o M i(s)-Ma(s) o 

V * * -(M2(*) + Gl)(iin-S)/ 

Point 3: for feedback law (11) /| can be written as follows, in (z,bi,b2,D) coordinates 

(-Hz{s) 0\ 

J* 2 =\ M i(S)-M2(i) 00 , 
\ * * M/ 

where M is a 2 x 2 matrix with the entries 

M= / -(M2(' f )+ G l)(' y in-'S) *in-« 

\-G 2 (fl2 (s) - D mia ) (£> max - M2 (s) ) 
Its eigenvalues are X\ = — ^(s) < 0, A2 = jUi (J) — jti2(*) < 0> A3 and A4 with 

A3A4 = G 2 (H2(s) - Anin)(Anax ~ /*2 (•*))( Jin ~ *) 

A3+A4 = -(At 2 (*) + Gi)(s i n-- f ) 

As in the proof of Proposition 2, one concludes the exponential stability of E\ when G2 > 
and Gi > — and the convergence of D(-) toward ^(s). □ 

Consequently, the adaptive control scheme proposed in Section 3 only allows one to 
reconstruct the larger of the two growth rates at any given s by stabilizing the equilibrium. 

Nevertheless, the introduction of feedback control may still be of help. To be specific, 
let us assume that species 1 is dominant for s < s c (and species 2 is suppressed there), 
and species 2 is dominant for s > s c , where s c is a cross-over point: fii(s) > [i 2 {s) for 
s < s c and [i\ (s) < [1 2 (*) for * > $c ( see the underlying function graphs in Figure 6 for a 
typical picture of the discussed scenario, s c = 0.5 in Fig. 6). Suppose we are interested in 
the location of E? to identify ^(s) for s < s c . As the form of the Jacobian J\ in (20) makes 
clear, two eigenvalues of /| are unaffected by our feedback control. One of them, — ^(f) is 
always stable. It corresponds to the transversally stable direction of the invariant subspace 
T = {(s,bi,b2) : s + b\ +b 2 = S{ n }. Once, the system is in T, control will not move it out of 
T, thus, we can ignore this eigenvalue. 

The other uncontrollable eigenvalue, [1\ (s) — ^(s), is unstable for s > s c , but stable for 
s < s c . Consequently, the following strategy would make it possible in principle to identify 



12 Jan Sieber et al. 



I±2 for s < s c : keep the system in the region where s > s c for some time to suppress species 
1 (which will exponentially decay for s > s c according to Proposition 3, point 2). When 
one is sufficiently close to the invariant line L2 = {(s,bi,b2) ■ b\ = 0,s + b2 = Sj n }, one is 
(nearly) in the single-species case, where one can then use the methods of Sections 3, 4 or 5 
to explore \Li for a finite time for s < s c until species 1 has recovered. This approach is also 
possible without control if both growth functions are monotone. 

Figure 6(b) demonstrates the application of feedback law (1 1) in combination with (14), 
which defines the adaptive control presented in Section 4, to the two-species situation with 
the monotonic Monod functions 

^W = frnr' ^ = L5 nLr (21) 

0.1+i (J.4 + J 

as growth rates. 

If one treats s as a parameter then the equilibria E\ and E\ undergo an exchange of 
stability (a degenerate transcritical bifurcation) at s = s c . As the adaptation rule (14) lets s 
drift slowly (with speed e) the full system exhibits a phenomenon known as delayed loss 
of stability [7,21] in the context of dynamic bifurcations [6], widely studied in slow-fast 
systems [22]. Physically this means (say, we are decreasing s slowly from below s c as in 
Figure 6(b)) that concentration b\ comes close to while s > s c . Then, when s has crossed 
s c , the concentration b\ grows exponentially, but still takes some time until it reaches values 
noticeably different from 0. The value si oss of fat which b\ becomes noticeably non-zero is 
in the ideal ODE model independent of the drift speed e of s. Figure 6(b) shows this effect: 
since b\ is nearly zero the variable D continues to follow the, by now unstable, drifting 
equilibrium £|. One can explicity determine the "input-output" map s h-> s\ oss for which one 
has s(t) ~ s and s(t + T) ~ $i oss at times t and t + T such that b\ (?) ~ b\ (t + T) : 

L? !oss 

when e is arbitrarily small. This delay mechanism allows one in principle a reconstruction 
of a part of the smaller growth rate close to the bifurcation value of s. 

Remark 1: From a practical view point, an apparent jump between the two graphs (see 
evolution curves in black in Fig. 6) could inform of the presence of another species, if one 
believes that the culture in chemostat was initially pure. Nevertheless, one can still rely on 
the reconstruction of parts of growth curves for each species. 




Fig. 6 Exploiting the delayed loss of stability for the dynamic control feedback law ( 1 1 ) in the two species 
case with growth rates given in (21) (shown underlying in both panels). The reference value f is increasing 
on the left panel, and decreasing on the right one. Parameters: e = 0.01 for panel (a) and e = —0.01 for panel 
(b) in (14), Gi = 2, G2 = 2 in (11), output disturbance defined in (13). 
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Fig. 7 Jumping back and forth between two regions at discrete times exploiting the delayed loss of stability 
for the simple control feedback law (2) in the two species case with growth rates given in (21) (shown 
underlying in both panels (c) and (d)). The reference value .v is increasing on the left panel, and decreasing 
on the right one. Parameters: Gi = 1, D milx = 1.5, D m in = 0.1 in (2), b m i n = 10~ 3 , output disturbance defined 
in (13). In panels (a-c) switching occurs between (s,D) = (0.1, 0.75) (where species 2 is stable) and (s,D) = 
(0.4 : 0.1 : 1, 0.9) (reading off jl^). In panels (d-f) switching occurs between (s,D) = (0.9,1) (where species 
1 is stable) and (s,D) = (1 : -0.1 : 0,0.15) (reading off /Ij). 
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Fig. 8 Effect of feedback control (2) on accuracy of identified growth rate jli . Parameters are the same as 
for Figure 7(a-c) (in particular, b mm = 10~ 3 ), the gain Gi is shown in the legend. The curves show the error 
relative to the difference between /ij and jli. For G\ = (no feedback control) D in (2) was varied to obtain 
different points approximating )![ , otherwise .v in (2) was varied (.v has no effect if Gj =0). 



Remark 2: In the idealized ODE model one could in principle recover the entire sup- 
pressed part of the growth rate by spending more time initially in the dominant part. This is 
so, because the concentration of the suppressed species (say, b\ for s > s c in Fig. 6(b)) can 
be made arbitrarily small by spending more time with s > s c . In practice, the suppression 
of species may not be perfect. For example, it may be impossible to suppress either species 
below a concentration fc m ; n > 0. Then this concentration b m [ n determines for how long the 
system will stay close to the invariant plane {bi = 0}, when this plane is unstable. Thus, Z? m j n 
determines how close one can get to the unstable equilibrium E£ for s < s c (together with the 
difference in growth rates, fi\ (s) — ^(.s), which determines how unstable the plane {b\ = 0} 
is in s). This is where the feedback control (2), D = sat^^ Dmax i (D — G\ (s — s)), has an ef- 
fect: the unstable equilibrium nas stronger attraction along the invariant line { (s , b 1 , bi ) : 
b\ = Q,s J rb\+b'i = ^n}, because the third eigenvalue of J\, — {^(s) + Gi )(■?;„ — f), can be 
made more negative by increasing G\ . 



14 



Jan Sieber et al. 



Figure 7 demonstrates that it is possible in principle to identify the growth rates of 
species even in ranges of s where they are suppressed, if the species is dominant in another 
region. The procedure was as follows (for Fig. 7(a-c)): 

1. Set (s,D) to (0.1,0. 75) , and wait until transients have settled (implying that species 2 
is suppressed). The output s settles to a value less than s c = 0.5. 

2. Then set (s,D) to, say, (0.9,0.75). One expects a transient that initially follows the 
invariant line {£>2 = 0,f>i = s ln — s] where s initially increases. 

3. As soon as 5 stops increasing (let's say, at .? out ), we know that the system now moves 
away from the plane {£>2 = 0}. So, we read off D, which is the estimate /Xi, e st(>Sout) (a 
black dot in Fig. 7(c)), and go back to step 1. 

Figure 7(d-f) demonstrates the same procedure for identifying for s < s c . The only dif- 
ference is that we read off D = ^estfsout) at an inflection point s oat of s(f). 

The procedure, with its steps 1-3 is also possible without feedback control (2). However, 
feedback control (2) increases the decay rate of the equilibrium with respect to disturbances 
within the invariant line {Z?2 = 0,bi = s m — s}. Thus, the system trajectory will come closer 
to the equilibrium before it diverges from the invariant line. Figure 8 demonstrates the effect 
of including the feedback term (2) if suppression of the unwanted species 2 is imperfect 
(^min = 10~ 3 ). The imperfect suppression is mimicked in our simulations by increasing £>, 
(i = 1,2) to 10~ 3 after each integration step if its value fell below 10~ 3 in this step. 

The two approaches in Fig. 6 and Fig. 7 correspond to two different choices for the 
trade-off between speed and accuracy. While the dynamic feedback in Fig. 6 requires only a 
single run, the procedure of switching back-and-forth between regions as rapidly as possible 
to obtain the growth rate of the suppressed species for abscissae s more distant from s c . 



7 Conclusion 

In this work, we have presented a new framework for the functional identification of non- 
monotonic growth functions in the chemostat. The proposed methods achieve identification 
by tracing out branches of equilibria also through their unstable parts. At the core of the 
method is the observation that the introduction of a stabilizing feedback loop transforms 
the problem of finding equilibria of the original uncontrolled (open-loop) system to a root- 
finding problem, which can then be solved using either continuous and step-wise variants 
of classical numerical continuation algorithms [2]. Numerical simulations illustrate the po- 
tential of the method. Further investigations are required to find to which extent it can be 
implemented in real experiments. 

The approach is more general than the case we have presented here for the chemostat 
model. We use the chemostat as a conceptually simple example that is still of practical 
interest. 

Another application we plan to explore in the future are regulation problems. For exam- 
ple, one can regulate the single-species chemostat to operate at the substrate concentration 
s at which the growth rate fi is maximal by following the same recipe. This approach to 
regulation, which is similar in spirit to the act-and-wait technique for delay compensation 
[17], does not require an a-priori identification of the growth rate /i, and leads to a different 
algorithm than the methods discussed in the literature [15, 16]). 
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